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ABSTRACT 

We study energy dissipation and heating by supersonic MHD turbulence in molecular clouds using 
Athena, a new higher-order Godunov code. We analyze the dependence of the saturation amplitude, 
energy dissipation characteristics, power spectra, sonic scaling, and indicators of intermittency in the 
turbulence on factors such as the magnetic field strength, driving scale, energy injection rate, and 
numerical resolution. While convergence in the energies is reached at moderate resolutions, we find 
that the power spectra require much higher resolutions that are difficult to obtain. In a 1024'^ hydro 
run, we find a power law relationship between the velocity dispersion and the spatial scale on which 
it is measured, while for an MHD run at the same resolution we find no such power law. The time- 
variability and temperature intermittency in the turbulence both show a dependence on the driving 
scale, indicating that numerically driving turbulence by an arbitrary mechanism may not allow a 
realistic representation of these properties. We also note similar features in the power spectrum of 
the compressive component of velocity for supersonic MHD turbulence as in the velocity spectrum of 
an initially-spherical MHD blast wave, implying that the power law form does not rule out shocks, 
rather than a turbulent cascade, playing a significant role in the regulation of energy transfer between 
spatial scales. 

Subject headings: ISM: clouds — ISM: magnetic fields — isothermal — simulations — stars: formation 
— turbulence 



1. INTRODUCTION 

Observed non-thermal line widths in molecular clouds 
(MCs), where all star formation in the Galaxy takes 
place, point to the presence of supersonic turbulence in 
such regions (Falgarone & Philips 1990). The proper- 
ties of the turbulent medium, such as Mach number and 
magnetic field strength, may determine the products of 
the star formation process. As an important source of 
heating within molecular clouds (Stone et al. 1998, here- 
after S98), turbulent energy dissipation may also play 
a role. Much effort has been directed toward numeri- 
cally simulating turbulent media in order to better un- 
derstand the link between turbulence and star formation 
(see Elmegreen & Scalo 2004, MacLow & Klessen 2004, 
and McKee & Ostriker 2007 and references therein). 

As we have not yet identified the turbulent driving 
mechanism, there remain many unanswered questions 
about the evolution of molecular clouds. Is the turbu- 
lence periodically re-energized, or does it simply decay 
away? How much impact do magnetic fields have on the 
properties of the turbulence? Crutcher (1999), using ob- 
servations of Zeeman splitting, found magnetic fields in 
some clouds strong enough that one cannot safely neglect 
their effects. Although it has been shown that magnetic 
fields do not appreciably lengthen the turbulent decay 
time (S98; Mac Low 1999), they do create anisotropy 
within the clouds (e.g. Vestuto et al. 2003, hereafter V03; 
Esquivel et al. 2003), which may have important obser- 
vational and evolutionary consequences. For example, 
molecular clouds are often observed to be filamentary 
(e.g. Mizuno et al. 1995; ChurchweU et al. 2004). 

In this paper, we will investigate the energy dissipa- 
tion properties of supersonic hydrodynamic and MHD 
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turbulence with a variety of magnetic field strengths us- 
ing Athena, a new higher-order Godunov code. An im- 
portant goal of this study is to investigate the effect of 
the assumed driving mechanism on the properties of the 
resulting turbulence, such as power spectra and intermit- 
tency indicators. Our analysis utilizes data from high- 
resolution numerical simulations with twenty-five differ- 
ent parameter sets. In recent years, a variety of results 
have been reported on the properties of supersonic MHD 
turbulence, including energetics (e.g. S98; Mac Low 1999; 
Ostriker et al. 2001), power spectra (e.g. Cho & Lazar- 
ian 2003, 2005; V03; Padoan et al. 2007, hereafter P07), 
and probability distribution functions (e.g. Padoan et 
al. 1997; Passot & Vazquez-Semadeni 1998; Ostriker et 
al. 2001; Kowal et al. 2007; Kritsuk et al. 2007, here- 
after K07). Where possible, we identify differences in 
our methods and results as compared to those of other 
groups. In a separate letter (Lemaster & Stone 2008, 
hereafter Paper I), we have reported the results of an 
investigation of the variation of the probability distribu- 
tion function (PDF) of density with Mach number. Our 
primary result in that paper was that the intermittent 
behavior of turbulence could be responsible for the large 
cloud-to-cloud variation in the observed star formation 
rate per solar mass. 

The primary tool available to investigate the proper- 
ties of highly nonlinear, supersonic MHD turbulence is 
direct numerical simulation. To date, most results have 
been computed using a few numerical algorithms, such as 
ZEUS (e.g. S98; V03; Mac Low 1999; Ostriker et al. 1999; 
Ostriker et al. 2001), the PENCIL code (e.g. Haugen 
& Brandenburg 2004), the Stagger code (e.g. P07), and 
ENO methods (e.g. Cho & Lazarian 2003). There has 
been concern expressed in the literature that previous 
results may be strongly affected by numerical dissipa- 
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tion. Thus, it is worth re-examining the problem with 
more accurate methods. This study represents one of 
the first apphcations of higher-order Godunov methods 
to the study of supersonic MHD turbulence. 

Without knowing the driving mechanism, we are left 
with many possible methods of generating turbulence in 
simulations. The hope is that the choice of method for 
the simulated driving will have an insignificant effect on 
the results. Federrath, Klessen, & Schmidt (in prep.), 
however, have shown that compressive and solenoidal 
forcing produce dramatically different turbulent states. 
For some diagnostics, such as intermittency, the time- 
variability of the turbulent state is critical. Even if an 
array of driving methods produce the same mean state, 
do the instantaneous states have the same distribution 
about the mean? We will investigate the dependence on 
driving scale of various diagnostics, given our particular 
driving method, which is very similar to that employed 
by, e.g., S98 and V03. 

For power spectra of simulated turbulence to be valu- 
able, the resolution needs to be high enough that the 
driving and numerical dissipation scales are well sepa- 
rated, allowing the inertial range to be studied. When 
magnetic fields are taken into account, simulating tur- 
bulence at these resolutions can be prohibitively expen- 
sive. Another goal of this study is to investigate whether 
properties of the power spectra and other diagnostics are 
converged at the numerical resolutions feasible at the mo- 
ment, up to 1024^. 

We begin in f|2] by describing our numerical methods 
in detail. We proceed, in fJ3J to present our results. H'3.l\ 
includes a convergence study and Mach number scaling 
study of saturation amplitudes, ^ 33. 21 includes a power 
spectrum analysis, §3.3l analvzes the sonic scale in 1024'^ 
runs, and !j3. 41 considers time- variability and temperature 
intermittency in the turbulence. Finally, we summarize 
our results in ij4] and discuss the implications. 

2. NUMERICAL METHODS 

The simulations we present were conducted up to a 
resolution of 1024^ with the Athena code (Gardiner & 
Stone 2005, 2008; Stone et al. 2008; Stone & Gardiner 
2008) on a three-dimensional Cartesian grid of length 
L = I with periodic boundary conditions. Athena uti- 
lizes a higher-order Godunov scheme which exactly con- 
serves mass, momentum, and magnetic flux. We solve 
the equations of ideal isothermal MHD, 

|f+V-(pv) = 0, (1) 

^-FV-(pvv-BB + P + BV2) = 0, (2) 
ot 

and 

— ^Vx(vxB), (3) 

where = 1 and P = c^p are the isothermal sound speed 
and pressure, respectively. Because our method of driv- 
ing the turbulence, described in the following paragraph, 
does not involve modifying the equations of MHD, we 
include no explicit forcing term here. We use an approx- 
imate nonlinear Riemann solver (HLLD; Miyoshi & Ku- 
sano 2005) for our MHD runs and an exact nonlinear Rie- 
mann solver for our hydrodynamic runs. Both our MHD 



and hydro simulations are integrated well past the tur- 
bulent saturation time using a directionally-unsplit van 
Leer scheme (Stone & Gardiner 2008). Further details 
of how we overcame the numerical difficulties associated 
with utilizing this method for turbulence can be found 
in Appendix [XI 

We initialize a uniform, stationary ambient medium 
with density p — 1 and magnetic field parallel to the 
X-axis whose amplitude Bq is fixed by the value of 
P = 2c^p/Bq. We then apply divergence-free velocity 
perturbations following a Gaussian random distribution 
with a Fourier power spectrum of the form 

|(5v^| cx fc6exp(-8fc/fcpk) (4) 

for kL/2TT < N/2, where N is the resolution and fcpk is 
the wavenumber of peak driving, in all but two of our 
runs. For the remaining two runs, with E/pL^c^ = 500 
and fcpki/27r = 2 at 1024^, we truncate the driving spec- 
trum at kL/2'K = 8. Before applying the perturbations 
to the grid, we shift them such that no net momentum 
will be added and normalize them to give the desired 
energy injection rate, E/pL^cl, or, in the decaying case, 
initial kinetic energy. 

For our driven runs, we choose our energy injection rate 
to give the desired turbulent Mach number, A4 = a^/cs, 
in the saturated state. Here, — [a^^ + a^^ + crlj^^^ 
is the 3D velocity dispersion of the gas and avi is 
the ID velocity dispersion in a given direction. The 
method we use in this paper to compute the turbu- 
lent Mach number differs from that used by some other 
groups; we utilize a density-weighted velocity dispersion, 

(^v, = [{pvf)/{p) - {{pv^)/{p))^]^/^. Note that, due to 
the zero net momentum of our turbulent medium, the 
second term on the RHS vanishes on the global scale. If 
we compare the Mach number computed with and with- 
out density-weighting the velocity dispersion, we find 
that the two are fairly correlated, but that the time- 
average of the latter is ~ 4% larger than that of the for- 
mer for driven strong- field MHD turbulence with, e.g., 
fcpkL/27r = 2 and E/pL'^cl = 500, or fcpki/27r = 8 and 
E/pL^cl — 1000. The comparison for the latter case is 
shown in Figure [TJ We will use a density- weighted veloc- 
ity dispersion for all but our sonic scale analysis. 

To drive turbulence in our simulations, we inject en- 
ergy before each time step, with a new realization of the 
power spectrum generated at intervals Atcs/L = 0.001. 
This differs from that done in S98 and V03, where the 
energy was injected only when a new realization of the 
spectrum was generated. Ostriker et al. (1999, 2001) 
and Kowal et al. (2007) used different driving spectra 
than our own but still constrained their velocity pertur- 
bations to be divergence free. Our method differs from 
that of K07 in that they approximate an isothermal equa- 
tion of state using an adiabatic index of 7 = 1.001, drive 
their turbulence on very large scales using an accelera- 
tion that is fixed in time, and also allow a substantial 
fraction of the energy introduced by the forcing to be in 
compressional modes. POT also used a fixed acceleration 
for turbulent driving. It is important to understand if 
these differences have any significant impact on the re- 
sults. 

For our decaying runs, we again begin with a uniform, 
stationary ambient medium but apply only one driving 
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Fig. 1. — Comparison of turbulent Mach number computed 
with (solid) and without (dashed) density- weighting for the driven, 
strong- field MIfD turbulence with k^]^L/2-K = 8 and E/pL^c'^ = 
1000. The mean of the latter (dash-dotted) is ^ 4% higher than 
that of the former (dotted). With the exception of for our sonic 
scale analysis, we use density- weighted velocity dispersions. 

impulse to the velocity field, with a power spectrum of 
the same form as is used to initialize our driven runs. 
We then allow it to evolve undisturbed. We choose the 
turbulent kinetic energy applied to our medium so as 
to give the desired initial turbulent Mach number. Our 
decaying runs differ from those of Sytine et al. (2000), 
who used a different initial driving spectrum and allowed 
a compressible component in the initial perturbations. 
They used an adiabatic equation of state with 7 = 1.4, 
which we will use only for our decaying simulations. In 
this case, we solve the total energy equation. 
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(S + P + BV2)v-B(B- v) =0, (5) 



in addition to Equations ([iJ-lS]), using the HLLC Rie- 
mann solver. Here, 
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(6) 



with P — (7 — l)e and = B • B, where e is the 
internal energy density and 7 is the ratio of specific heats. 
The numerical methods in Athena conserve total energy 
exactly. 

We primarily investigate strong-field MHD turbulence 
with (3 = 0.02 and hydrodynamic turbulence (/? = 00), 
although we present runs with (3 = 0.2 (moderate- 
field) and (3 — 2.0 (weak-field) as well. Note that, 
due to a definition of (3 in S98 and V03 that differs 
from ours by a factor of 2, our (3 — 0.02 results should 
be compared to their j3 = 0.01 results, and similarly 
for other values of (3. The magnetic fields we use in 
our simulations correspond to physical values of i? = 
2.0//G/3-i/2(r/10K)i/2(^^jl02cm-3)i/2^ where T is 
the temperature and n//^ is the number density of molec- 
ular hydrogen. Our simulations are scale-free, allowing 
them to be scaled to any set of physical parameters using 
appropriate choices of p, Cs, and L. Utilizing the same 
values given in S98, i.e. L = 2pc, n/f^ = lO'^cm"^, 



and T ~ 10 K, yields energy injection rates of up to 
E — 0.4^0, with a magnetic field strength B = 44 /iG 
for the strong-field case. 

3. RESULTS 

Figure [5] shows slices in mass density along the far 
faces of the computational domain for driven = 6.9 
strong-field MHD turbulence with fcpkL/27r = 2 at 1024^. 
Also included are magnetic field vectors in a slice normal 
to the y-axis near the bottom of the cube. Due to the 
strong background magnetic field, the vectors are fairly 
well aligned. Figure [3] also shows slices in density, but 
this time for M — 7.0 hydrodynamic turbulence with 
the same driving scale and resolution. For both cases 
the structure seen is filamentary, with anisotropy in the 
scale of the structures in the strong-field case as a result 
of the field. This is not to say, however, that the fila- 
ments are aligned with the magnetic field. In fact, some 
appear to be oriented perpendicular to the field. Figures 
2] and \5\ show column density along the line of sight par- 
allel to the z-axis for the same runs as in Figures [2] and 
[31 respectively. The projected filamentary structure is 
visible in these column density images. The column den- 
sity contrast is higher for the strong-field than for the 
hydrodynamic case. 

3.1. Saturation Amplitude 

We begin our quantitative analysis by studying the en- 
ergy in fluctuations once our driven turbulence runs have 
reached saturation. Since the method with which we 
drive our turbulence injects energy at a constant rate, 
at saturation the energy dissipation rate of the turbu- 
lence will, on average, equal the energy injection rate. 
At sufficiently high numerical resolution, the numerical 
dissipation will become negligible compared to the phys- 
ically interesting sources of dissipation (shocks), and the 
energy dissipation properties of the turbulence can given 
us insight into the heating within and evolution of molec- 
ular clouds. In this section, we will first determine the 
resolution at which the saturation energy has converged, 
i.e. reached a state such that further increasing the reso- 
lution has a negligible effect on the state, for turbulence 
evolved using Athena and compare the turbulent decay 
rates to those presented in S98 computed using ZEUS. 

The kinetic energy associated with the fluctuations 
in our turbulent medium can be quantified by Ek = 
Jpw^/2(ix, where p is the density and v is the mag- 
nitude of the velocity. An integration for this and all 
similar quantities is performed over the entire domain. 
Similarly, the energy in magnetic field perturbations is 
SEb = J{B^ — -Bo)/2dx, where B and Bq are the mag- 
nitudes, respectively, of the magnetic field and its ini- 
tial value. The total energy in perturbations is the sum 
of the energy in kinetic and magnetic field fluctuations, 
Pport = Ek + SEb- To analyze the energy dissipa- 
tion properties of the turbulence, we compute dissipation 
timescales using idiss — E/E, normalizing them to the 
flow crossing time, t{ ~ Apk/ \/2Ek, where Apk = 27r/fcpk 
is the wavelength of peak driving. 

We can partition the kinetic energy by breaking down 
the velocity into its solenoidal and compressive com- 
ponents. The solenoidal component is divergence- free, 
whereas the compressive component is curl-free. These 
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Fig. 2. — Driven strong-field M ~ 6.9 MHD turbulence with = 2 at 1024'^. Slices of density along the far faces of the cube 

on a logarithmic color scale from 0.02 (blue) to 5.0 (red). Magnetic field vectors along a slice normal to the y-axis aX, y ^ 0.0625 are fairly 
well aligned and there is anisotropy in the scale of the structures that results from the magnetic field. 

to 512'^. We study two sets of simulations with identi- 
cal energy injection rate, E/pL^cl = 1000, but differing 
driving scales. The set of runs with small-scale driving, 
fcpki/27r — 8, correspond to the hydrodynamic run in 
S98, while the other set of runs are driven at twice the 
scale, fcpki/27r = 4. The properties of the 512^ run from 
each set can be found in Table [TJ We have not performed 
a convergence study of the fcpki/27r = 2 runs due to the 
high level of time- variability (discussed in ! j3.4.ip . 

Although both sets of runs are driven with energy in- 
jection rate E/pL'^c^ = 1000, the larger driving scale of 
the fcpki/27r = 4 set causes it to converge to « 7.2, 
while the fcpki/27r = 8 set converges to only Ai ~ 5.8. 
The former set reaches convergence by 64^, with higher 
resolutions having a small scatter about the converged 
value. The latter set, on the other hand, converges mono- 
tonically at order a ~ 1.6 for all resolutions analyzed, 
reaching within 1% of the converged value by 128'^. 

The kinetic energy in fluctuations of the set with larger 
driving scale, kp^L /2Tr = 4, converges to a value of Ek « 
26 by 64^^ and has some scatter about the converged value 
for higher resolutions. For the set with smaller driving 
scale, fcpk£/27r = 8, the kinetic energy in fluctuations 
converges at order a « 1.5 to only Ek ~ 17. The value 
at 128^ is within 2% of the converged value, while at 256^ 



can easily be computed using V(7(k) = [k ■ v(k)]fc and 
V5(k) = [fc X v(k)] X k, respectively. For simulations 
with comparable parameters (i.e. kp\^L/2Tr = 8 and 
E/pL^cl = 1000), these energies can be compared di- 
rectly with the values given in S98 and V03 after account- 
ing for the difference in definition of /3. We average our 
quantities over at least one dynamical time, often several, 
beginning after the turbulence has fully saturated. 

To determine the rate of convergence of a given volume- 
integrated quantity, we compute the percent error in 
the quantity at each resolution relative to the converged 
value. In the case where the quantity changes mono- 
tonically with resolution, we find the converged value by 
performing a three-parameter fit, [qn — <loo) / <lcx, — cN^", 
where N is the resolution and g^v is the value at that res- 
olution. The results of this fit tell us (1) the converged 
value, goo, (2) the order of convergence, a, and (3) the 
coefficient, c, that determines the resolution at which our 
result has converged, i.e. when (qjv — '?oo)/<Zoo < 0.01. 

3.1.1. Hydrodynamic Convergence Study 

To determine the resolution at which our numerical 
dissipation has become small compared to shock dis- 
sipation, we analyze the properties of driven, super- 
sonic hydrodynamic turbulence at resolutions from 32^ 
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TABLE 1 
Driven Hydro Turbulence at 512'^ 







M 


Ek 


Ec 


Es 


idiss/tf 


ge/E 


<^Q,v/Qv 


o'q,m/Qm 


8 


1000 


5.8 


17 


3.8 


13 


0.77 


< 1% 






4 


1000 


7.2 


26 


5.9 


20 


0.76 


1% 


1% 


2% 


4 


375 


5.3 


14 


3.1 


11 


0.78 


1% 


1% 


2% 


4 


140 


3.8 


7.2 


1.6 


5.6 


0.78 


1% 


2% 


2% 


4 


40 


2.6 


3.3 


0.7 


2.6 


0.84 


1% 


< 1% 


< 1% 


4 


3.5 


1.2 


0.8 


0.1 


0.7 


1.1 


< 1% 


< 1% 


< 1% 


2 


500 


7.0 


25 


5.5 


19 


0.7 


4% 


2% 


4% 


2 


187.5 


5.2 


13 


3.1 


10 


0.73 


2% 


1% 


4% 


2 


70 


3.7 


6.9 


1.6 


5.4 


0.74 


2% 


2% 


2% 


2 


20 


2.5 


3.1 


0.6 


2.5 


0.77 


1% 


3% 


3% 


2 


1.75 


1.2 


0.7 


0.1 


0.6 


0.99 


2% 


4% 


4% 



it is within a fraction of a percent. 

Although we have an extremely low sampling rate for 
the fraction of the kinetic energy in solenoidal and com- 
pressive modes, they appear to be independent of driving 
scale. The solenoidal fraction converges to Es/Ek ~ 
0.78, decreasing with increasing resolution, while the 
compressive fraction converges to Ec/Ek ~ 0.22. The 
solenoidal fraction is within a fraction of a percent of 
the converged value by 256'^; however, it varies by only 
a small amount down to low resolution. Finally, the ra- 
tio of the energy dissipation timescale to the flow cross- 



ing time at the driving scale increases with resolution, to 
idiss/if ~ 0.78 for the fcpki/27r = 4 set and tdiss/tf ~ 0.76 
for the kpyL/2Ti — 8 set, only a small difference. Con- 
sistent with our previous results, the k^^^Ljlix — 4 set 
converges at lower resolution than the k-p]^L/2TT = 8 set, 
64'^ and 256'^, respectively. The turbulence dissipates in 
less than a flow crossing time in all cases. 

Although the converged Mach number and saturation 
energies are higher for fcpki/27r = 4 than for k-pyL/2'K = 
8, the fraction of the energy in the solenoidal or compres- 
sive mode, as well as the ratio of the energy dissipation 
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Fig. 4. — Column density along line of sight parallel to z-axis for driven strong-field MHD turbulence with kp]^L/'2n = 2 at 1024'^ on a 
linear gray scale from 0.3 (black) to 3.4 (white). 

timescale to the flow crossing time at the driving scale, 
are relatively independent of the driving scale. While 
most of the quantities of interest converge between 128'^ 
and 256'^ for the kp]^L/2Tr = 8 hydrodynamic runs, con- 
vergence has already been reached for these quantities by 
64"^ for the k-p\^L/2'K ~ 4 runs. High resolutions are crit- 
ical for separating the driving and dissipation scales in 
the power spectra far enough to study the inertial range; 
however, quantities such as energy dissipation rate and 
turbulent Mach number converge at resolutions which 
are more easily attainable. 

3.1.2. MHD Convergence Study 

We now analyze driven, strong-field, supersonic MHD 
turbulence in the same manner as in the previous section. 
To determine the resolution at which our numerical dis- 
sipation becomes negligible, we study the convergence of 
two sets of MHD simulations with energy injection rate 
E/pL'^c?^ = 1000 at resolutions from 32^ to 5123. Qne set 
has k-pyL/2TT = 8, similar to the strong-field MHD run in 



S98, while the other has kp\^L/2-K — 4. The properties 
of the 512'^ run from each set can be found in Tabled 
Again, we have not performed a convergence study of the 
fcpki /27r = 2 runs due to the high level of time- variability 
(discussed in ii3.4.ip . 

As in the hydrodynamic case, the fcpki/27r — 4 MHD 
runs converge to a higher Mach number, Ai ~ 6.8, than 
the fcpki/27r = 8 runs ~ 5.5). The former is con- 
verged by 256^, while the latter converges at I*'' order, 
only reaching convergence at 512^. The total energy in 
fluctuations converges at nearly P* order to -Epcrt ~ 36 
for fcpkL/27r = 4, shown in Figure El While convergence 
has definitely been reached by 512'^, the 256^ value is less 
than 2% from the converged value. For kp\iL/2Tr = 8, on 
the other hand, convergence is not reached until 512'^, 
where import ~ 23. Convergence is at nearly I*'* order for 
this case as well. 

For the energy in magnetic field perturbations, con- 
vergence is reached at 512^ for both the fcpki/27r — 4 
and fcpki/27r = 8 cases. Convergence is approached at 
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Fig. 5. — Column density along line of sight parallel to z-axis for driven hydrodynamic turbulence with kp\^L/2'K = 2 at 1024^ on a 
linear gray scale from 0.3 (black) to 3.4 (white). 



TABLE 2 
Driven MHD Turbulence at 512^ 



p 


fcpkL/27r 




M 


Ek 


Ec 


Es 


Eb 


^diss/^f 


<Je/E 




<yQ,M IQm 


2.0 


8 


1000 


5.3 


14 


2.1 


12 


5.8 


0.83 


< 1% 






0.2 


8 


1000 


5.1 


13 


1.6 


11 


8.6 


0.89 


< 1% 






0.02 


8 


1000 


5.4 


15 


1.5 


i:-! 


8.2 


1.0 


< 1% 






0.02 


4 


1000 


6.8 


23 


2.3 


21 


13 


0.98 


< 1% 


1% 


2% 


0.02 


4 


375 


5.2 


14 


1.1 


13 


6.4 


1.1 


< 1% 


< 1% 


2% 


0.02 


4 


140 


3.8 


7.3 


0.5 


6.7 


3.1 


1.1 


3% 


2% 


< 1% 


0.02 


4 


40 


2.7 


3.6 


0.2 


3.4 


1.4 


1.3 


< 1% 


3% 


3% 


0.02 


4 


3.5 


1.3 


0.8 


0.04 


0.8 


0.3 


1.7 


< 1% 


2% 


2% 


0.02 


2 


500 


7.0 


24 


2.2 


22 


11 


0.98 


4% 


4% 


5% 


0.02 


2 


187.5 


5.0 


13 


1.2 


11 


5.9 


0.99 


3% 


2% 


2% 


0.02 


2 


70 


3.7 


6.8 


0.6 


6.3 


2.9 


1.0 


3% 


4% 


5% 


0.02 


2 


20 


2.5 


3.1 


0.2 


2.9 


1.3 


1.1 


2% 


3% 


4% 


0.02 


2 


1.75 


1.2 


0.8 


0.06 


0.7 


0.4 
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Fig. 6. — Total energy in fluctuations in driven MHD turbulence 
runs with = 4 at 32^ (dotted), 64^ (short dashed), 128^ 

(long dashed), 256^ (dash dotted), and 512'^ (solid). The saturated 
energy has converged by 256'^. 

roughly 0.7 order for both cases. The fraction of the total 
energy (magnetic plus kinetic) in magnetic field fluctu- 
ations increases with resolution, reaching 35% at 512'^. 
When driven at fcpki/27r — 4, the kinetic energy in fluc- 
tuations converges to Ek ~ 23 by around 256'^, while at 
fcpki/27r = 8 it converges at order 0.9 to Ek ~ 15 by 
5123. 

The fraction of the kinetic energy fluctuations in 
solenoidal modes increases with resolution for both the 
fcpki/27r — 4 and kp]^L /2tt ~ 8 cases, converging to 
Es/Ek ~ 0.9. The fraction in compressive modes, on 
the other hand, decreases with resolution to Ec/Ek ~ 
0.1. Convergence for these quantities is not reached un- 
til 512'^. Finally, the ratio of the dissipation timescale 
to the flow crossing time at the driving scale converges 
to tdiss/if ~ 0.98 by 256^ for fcpkL/27r = 4, while for 
fcpki'/27r = 8 it converges to tdiss/if ~ 1-0 by 512"^. Even 
with a strong magnetic field, the turbulence dissipates in 
a fiow crossing time. 

Just as in the hydrodynamic case, the increased driv- 
ing scale caused a higher turbulent Mach number and 
energy at saturation for the fcpki/27r = 4 than for the 
fcpki/27r = 8 case. The fraction of the kinetic energy 
in solenoidal or compressive modes, however, was inde- 
pendent of driving scale, just as before. The ratio of 
the dissipation timescale to fiow crossing time was also 
relatively independent of the driving scale. Convergence 
for strong-field MHD turbulence was reached at higher 
resolutions, typically 512'^, than for hydrodynamic tur- 
bulence. 

3.1.3. Comparison to ZEUS 

Because we performed a convergence study using the 
same turbulence parameters as in S98, we can directly 
compare our high-resolution results to those published 
therein. For driven strong-field MHD turbulence with 
fcpki/27r = 8, we find at 256^ that our total energy in 
fluctuations at saturation is only 8% higher than that 
found in S98, due to the lower level of numerical dissi- 
pation at this resolution in Athena than in ZEUS. For 



intermediate- and weak- field MHD, our energies are 7% 
and 8% higher, respectively. It is likely that neither our 
results nor the S98 results have converged by this reso- 
lution, however. 

At 512^, there is no obvious difference between the to- 
tal energy in fiuctuations for our strong-field MHD sim- 
ulation and the S98 result. The values from these two 
codes converge, even though they utilize completely dif- 
ferent numerical methods. The ratio of magnetic to ki- 
netic energy fluctuations, however, is different for the 
two codes, yielding ratios of dissipation timescale to flow 
crossing time that differ more substantially. At 256'^, this 
ratio is 14% greater for our strong-field MHD run than 
that presented in S98. For intermediate- and weak-field 
MHD, our ratios are both 11% greater. In all cases, the 
ratio of timescales remains below unity. 

For our driven hydro run at 256^ with kp]^L/2'K = 8, 
we find that our total energy in fluctuations at saturation 
is also only 8% higher than that found in S98. Although 
our result has converged by this resolution, the ZEUS 
result may not yet be converged. The ratio of dissipa- 
tion timescale to fiow crossing time in this case is 11% 
greater than that found in S98, but is still below unity. 
Although it has been suggested in the literature that the 
rapid decay of supersonic turbulence is due to excessive 
numerical dissipation in ZEUS, these results, computed 
with a higher-order Godunov scheme, do not support 
that conclusion. 

3.1.4. Hydrodynamic Mach Number Scaling 

To investigate the effect of the turbulent Mach num- 
ber on energies and dissipation rates in the turbulence, 
we now analyze two series of five driven, supersonic hy- 
drodynamic turbulence simulations at 512'^. One series 
is driven at fcpki/27r = 4, while the other is driven at a 
larger scale of fcpki/27r — 2. The energy injection rates 
of the latter series, the same that we use in Paper I, 
are half that of the former series, giving roughly equal 
Mach numbers to the corresponding pairs. These Mach 
numbers all fall within the range 1.2 < Ai < 7.2. The 
properties of these runs can be found in Table [TJ 

The kinetic energy in each series of runs spans 1.5 or- 
ders of magnitude. We find a power law relationship 
between the kinetic energy and energy injection rate, i.e. 

Ek « O.A9[ikp^L/27:)/2]-'^^E/pL'cl)°-'\ (7) 

The relationship between the kinetic energy and Mach 
number is, of course, exactly Ek = 0.5A^^ since our 
Mach numbers are computed from the density-weighted 
velocity dispersion and there is no net momentum asso- 
ciated with the turbulent medium as a whole. We find 
for our driving method that we can estimate the Mach 
number that will result from a given energy injection rate 
using 

M » Q.99[{kpi,L/2TT)/2]-^/\E/pL\lf-^\ (8) 

While this equation gives an indication of how the dissi- 
pation scales with Mach number, the exact relations will 
only apply to our unique driving method. 

For the runs driven at kp]^L/2'K — 4, the fraction of 
the kinetic energy in compressive modes varies by only a 
fraction of a percent among the runs with 3.8 < < 7.2. 
These runs have Ec/Ek ~ 0.22. The fraction drops 
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off only slightly for the A4 — 2.6 run, but then drops 
substantially, to Ec/Ek ~ 0.14, for the M = 1.2 run. 
For the runs driven at fcpk-L/27r = 2, the fraction for the 
M = 7.0 and M = 3.7 runs is also Ec/Ek « 0.22, while 
the M — 5.2 and Ai —2.5 runs have a bit larger and 
smaller fractions, respectively. Once again, the fraction 
for the = 1.2 run is substantially smaller, Eq/Ek ~ 
0.15. It appears as though the fraction of the energy 
in compressive modes remains roughly the same except 
when the turbulence is only mildly supersonic, where the 
fraction is much lower. 

The ratio of the dissipation timescale to the flow cross- 
ing time at the driving scale remains below unity for all 
except the M. = 1.2 run at kp]^L/2TT = 4. While the ratio 
for the runs at kpkL/2'K = 4 with M > 3.8 doesn't vary 
much, when we increase the driving scale to kp]^L/2'K = 2 
the value does change with Mach number. The ratio for 
the k.pyL/2TT = 4 runs is always larger for a given Mach 
number than is the ratio for the kpi^L/2T: = 2 runs. 

3.1.5. MHD Mach Number Scaling 

Finally, we analyze two series of five driven, supersonic, 
strong-field MHD turbulence simulations at 512'^ to de- 
termine the effect of the turbulent Mach number on the 
energies and dissipation rates. As before, one series is 
the same kp]^L /2Tr = 2 set of runs used in Paper I, while 
the other is driven at k-p\^L/2'K = 4 with twice the energy 
injection rate, yielding pairs of simulations of roughly 
the same Mach number within the range 1.2 < M < 7.0. 
The properties of these runs can be found in Table [H 

As in the hydro case, the kinetic energy in each strong- 
field MHD series spans nearly 1.5 orders of magnitude. 
The total energy in fluctuations, which also includes the 
energy in magnetic field perturbations, increases by an 
even larger amount from the lowest to highest Mach num- 
ber runs. The power law relationship between the Mach 
number and energy injection rate is 

M « (1.04±0.02)[(fcpki/2^)/2]-i/4(ij/pL2c3)0-30, (9) 

which leads to a slightly narrower range of Mach numbers 
than does the relation for pure hydro. Again, although 
this gives an indication of how energy dissipation scales 
with Mach number, the exact relationship is unique to 
our driving method. 

The fraction of the kinetic energy in compressive modes 
for strong-field MHD is much less than for hydro. When 
fcpki/27r — 2, the fractions increase monotonically with 
Mach number, from only Ec/Ek ~ 0.05 for M — 1.3 
to Ec/Ek « 0.10 for M = 6.8. When fcpkL/27r = 4, 
however, the values are much less predictable. In this 
case, the = 3.7 and = 1.2 runs have nearly iden- 
tical values when kp\^L/2'K = 2, whereas for hydro the 
fraction at low Mach number was substantially smaller 
than the values for higher Mach numbers. 

The ratio of the dissipation timescale to the fiow cross- 
ing time at the driving scale decreases with Mach num- 
ber. While the values for high Mach number are below 
unity, the low Mach number values become as large as 
tdiss/U ~ 1-7. The values for the highest Mach numbers 
are roughly the same for both driving scales, but for the 
smaller driving scale they increase more quickly. The 
Mach number has a strong infiuence on the dissipation 
rate of the turbulence, but even for w 1.2 the dissi- 



pation timescale does not exceed twice the flow crossing 
time. 

3.2. Power Spectra 

We next consider turbulent velocity power spectra. 
We compute the power spectrum (PS) of the veloc- 
ity, -PK(k) — |v(k)p/2, as well as that of its compres- 
sive and solenoidal components, -Pc(k) — \k ■ v(k)p/2 
and -Ps(k) = |A; x v(k)p/2, respectively, in the same 
way as was done in V03 for the PSD of the specific 
kinetic energy and its components. To generate the 
spherically-integrated compensated power spectra that 
we will present, we average /'(k) within spherical shells 
and multiply by the volume within the shell, dV = 
Airk^dk, to find P{k), where k ^ {k^ + k^ + klfl"^. 

To determine the degree of anisotropy in MHD turbu- 
lence, we will also analyze cylindrically-integrated power 
spectra. We generate these spectra by averaging Pk{^) 
over cylindrical shells whose axes are oriented along the 
mean magnetic field direction. This yields , 
where fcy = \kx\ is parallel to the mean magnetic field 
and fc_L = {ky -\- k^Y^"^ is in the plane normal to the field. 

The majority of the power spectra we present in our 
figures are compensated — divided by a power law to pro- 
duce plots where the inertial range is very roughly hori- 
zontal, making small deviations from a power law easier 
to see. For each of our driven turbulence runs, we com- 
pute power spectra from many snapshots taken at reg- 
ular intervals, averaging them together to minimize the 
effects of rare events. With the exception of the 1024'^ 
runs, we average over at least 69 snapshots spanning at 
least 4.6 dynamical times. While we averaged over 2.7 
dynamical times for the 1024^^ hydro run, the MHD case 
is not sufficiently time-averaged. 

3.2.1. Decaying Subsonic Hydrodynamic Turbulence 

Sytine et al. (2000) presented the power spectra of de- 
caying, subsonic, adiabatic hydrodynamic turbulence at 
a range of resolutions. These power spectra demonstrate 
the formation of a feature known as the bottleneck ef- 
fect. Energy cascades down from larger scales but can- 
not easily be dissipated, causing a build-up of small-scale 
power. We begin by verifying that Athena can reproduce 
the bottleneck. These are the only set of adiabatic runs 
we will consider. 

The initial velocity perturbations in these decaying 
runs have the same form of power spectrum as in our 
driven runs with kp\^L/2TT = 4 and are normalized to 
give an initial turbulent Mach number of Mo = 0.5. 
After this impulse is given to the initially uniform am- 
bient medium, it is allowed to evolve undisturbed until 
tCs/L = 2. Figure [7] shows the uncompensated veloc- 
ity power spectrum at this time for runs with resolu- 
tions from 256^^ to 1024'^. As expected, this plot ap- 
pears very similar to Figure 11 of Sytine et al. (2000). 
Since Athena was designed to have low numerical dis- 
sipation, the bottleneck is strong in our simulations of 
mildly-compressible (subsonic) turbulence. 

Analyzing the compressive and solenoidal components 
of the velocity separately, we find that the bottleneck 
appears quite strong in the PS of the latter while appar- 
ently absent in that of the former. This is not surprising 
as shocks directly couple large and small scales, allowing 
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Fig. 7. — Spherically-integrated velocity power spectra at 
tCs/L = 2 for decaying initially Mach 0.5 adiabatic hydro tur- 
bulence with kp^L/2iT = 4, at 1024^ (solid), 512^ (long dashed), 
and 256^ (dotted). Compare to Figure 11 of Sytine et al. (2000). 

the compressive energy to bypass the turbulent cascade 
and be immediately dissipated. The absence of the bot- 
tleneck in the compressive component agrees with the re- 
sults of Porter et al. (1999) for Mach 1 driven, adiabatic 
hydro turbulence. At the very highest wavenumbers in 
our simulations, the power in the compressive compo- 
nent flattens out. Although this could be an artifact of 
the code related to its treatment of shocks, one should 
be wary of putting too much stock in the high-/c region 
of any power spectrum as effects such as aliasing (in- 
troduced in the calculation of the power spectrum, not 
the fluid dynamics) could substantially alter the power in 
that region (see, e.g., Press et al. 1992 for a more detailed 
discussion) . 

3.2.2. Driven Hydrodynamic Turbulence 

Next we consider driven supersonic (isothermal) hydro- 
dynamic turbulence. We compute velocity power spectra 
for runs with k-p\^L/2'K = 2 to maximize the separation 
between the driving and dissipative scales. To obtain a 
slightly higher turbulent Mach number, M « 7.0, than 
before, we use an energy injection rate of E/pL^C^ — 500. 
These are the same runs from the Mach number scaling 
analysis of ^13.1.41 

Figure [5] shows our compensated time-averaged veloc- 
ity power spectrum for resolutions of 256^ through 1024^. 
To align the dissipative range for simulations with dif- 
ferent resolutions, we express wavenumber as a fraction 
of the Nyquist value, k^^L /2Tr = N/2 for a simulation 
with resolution N^. At low resolution we see no iner- 
tial range, but by 1024^^ we have separated the driving 
and dissipative scales enough that we may be just start- 
ing to see an inertial range. There is a small range, 
6 < kL/2TT < 11, where we see roughly a P{k) oc fc"^ 
power law, but we are still using a small amount of forc- 
ing up through kL/2'K = 8 so the slope in this range may 
be affected. We find a much longer power law range over 
11 < kL/2TT < 36, P{k) cx fc^i '^. Although it may be 
argued that a slope this shallow must be due to a bot- 
tleneck, the range of scales with steeper slope is far too 
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Fig. 8. — Spherically-integrated, compensated power spectra 
of velocity for driven hydro turbulence runs with fcp]jL/27r = 2 at 
1024^ (solid), 5123 (jo^g dashed), and 256^ (short dashed). The 
X-axis has been renormalized to give fc/fcjv = 1- Also shown are fits 
(dotted) to the slope of the 10243 run over the ranges 6 < kL /2n < 
11 and 11 < kL/2iT < 36, representing power laws P{k) oc k~^-'^^ 
and P{k) oc k~^'^'^, respectively. 

limited to conclude this definitively from our data. A 
higher resolution would be required to determine if this 
is actually the case. 

Figures and [10] show power spectra of the solenoidal 
and compressive components of velocity, respectively. 
For the solenoidal component, we find power laws P{k) oc 
and P{k) cx k-'^-^ over the ranges 6 < kL/2'n- < 11 
and 11 < kL/2-K < 36, respectively. For the compressive 
component, however, we find only P{k) oc fc~^ ® over 
the range 10 < kL/2'K < 40. In two independent 512'^ 
simulations that use different random perturbations to 
seed and drive the turbulence, we find the slopes of the 
compressive component to agree with each other as well 
as that of the 1024'^ simulation to within a percent, sug- 
gesting that the length of our time- averaging is sufficient. 
The ratio of power in the compressive and solenoidal 
components of velocity is much higher than for subsonic 
turbulence, having the effect of slightly washing out the 
shallow feature in the velocity spectrum in the supersonic 
case. 

When comparing hydro turbulence with kp]^L/2n = 8, 
we find more power at high wavenumbers for turbulence 
evolved using Athena than using ZEUS. Although it was 
stated in V03 that the bottleneck would not affect fi- 
nite difference codes, plotting the compensated power 
spectrum shows some evidence of blending of the driving 
peak with a secondary bump or shallowing of the spec- 
trum (much less prominent than in our own). If this is 
indeed a bottleneck, this difference likely results from the 
greater numerical dissipation in ZEUS than in Athena. 

3.2.3. Driven MHD Turbulence 

Analyzing driven strong-field MHD turbulence with 
fcpki/27r = 2 and an energy injection rate of E/pL'^c^ — 
500, shown in Figure fTTl we again find no inertial range at 
low resolution. In this Ai ~ 7.0 run, the velocity power 
spectrum is dominated by the solenoidal component, re- 
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Fig. 9. — Spherically-integrated, compensated power spectra of 
the solenoidal component of velocity for driven hydro turbulence 
runs with = 2 at 1024^ (solid), 512^ (long dashed), and 

256^ (short dashed). The X-axis has been renormalized to give 
k/kff = 1. Also shown are fits (dotted) to the slope of the 1024'^ 
run over the ranges 6 < kL/2-K < 11 and 11 < fcL/27r < 36, repre- 
senting power laws P{k) oc fc"^"'* and P{k) oc k~^'^^, respectively. 
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Fig. 10. — Spherically-integrated, compensated power spectra of 
the compressive component of velocity for driven hydro turbulence 
runs with fcp]jL/27r = 2 at 1024^ (solid), 512"^ (long dashed), and 
256^ (short dashed). The X-axis has been renormalized to give 
k/kff = 1. Also shown is a fit (dotted) to the slope of the 1024'' 
run over the range 10 < kL/2n < 40, representing the power law 
P{k) oc fc-l-Sl. 

suiting in a substantial shallowing of the spectrum that 
is apparent in the higher-resolution runs. At 512"^, the 
slope of the spectrum over the interval 8 < kL/2-K < 18 is 
only slightly steeper than P{k) cx k~*^^ for the velocity 
and its solenoidal (Figure fT2|) component. We find spec- 
tral slopes in two independent 512"^ simulations (whose 
power spectra were averaged together to give that plot- 
ted in the figures) that agree to within a few percent, 
suggesting that our time-averaging to obtain a reliable 
spectral slope. 
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Fig. 11. — Spherically-integrated, compensated power spectra 
of the velocity for driven strong-field MHD turbulence runs with 
kpi,L/2n = 2 at 1024^ (solid), 512^ (long dashed), and 256^ (short 
dashed). The x-axis has been renormalized to give fc/fcjv = 1. Also 
shown is a fit (dotted) to the slope of the 512'^ run over the range 
8 < fcL/27r < 18, P{k) oc The 1024^ run appears only 

slightly steeper. 

Although we do not have a sufficiently long time- 
average to find a robust fit to the slope of the 1024^^ 
run, the slopes appear to be very similar. As before, 
the slope we have found is much shallower than is typi- 
cally predicted predicted for the inertial range, while be- 
ing steeper than that presented in K07 as the bottleneck 
spectral slope. We note that in the MHD case, the range 
over which we find a power law seems to increase with 
resolution, which one would not expect if the feature was 
due to numerical dissipation (bottleneck). 

It has been argued by K07 that the quantity u = p^/'^v 
should have an inertial range power law of k~^^^ even 
for supersonic turbulence, but we see over a wide range 
of wavenumbers a k~^/^ law (Figure [T5|) . Kurien et 
al. (2004) have shown for subsonic turbulence that, while 
the k~^/^ law is appropriate when the energy timescale 
dominates over the helicity timescale, a law is ex- 

pected when the helicity timescale is non-negligible com- 
pared to the energy timescale, even when the relative he- 
licity may seem insignificant. This should be investigated 
as a possible cause of the bottleneck effect in supersonic 
turbulence. 

Over the range 10 < kL/2'K < 18, the compressive 
component (Figure fT4|) of the 512^ run follows a power 
law P{k) oc k~^'^. For 1024^^ run, the power spectrum 
appears to smoothly change slope, indicating that it is 
not a power law, but a longer time-average might change 
the shape. The power spectrum of the magnetic field at 
512'^ follows the power law P{k) cx k~^-^ over the range 
8 < kL/2'K < 18. At 1024"^, the slope appears very simi- 
lar, but with extra noise at the smaller wavenumbers due 
to the short time-average used. For all resolutions stud- 
ied, the uncompensated power approaches a constant, 
non-zero value at the highest wavenumbers. This pro- 
duces a very prominent upturn in the compensated spec- 
tra. 

A direct comparison for strong-field MHD with 



12 



Lemaster & Stone 



10-2 







1 1 




: 










/ / 


/ 




- 


/ / 


/ 


\\\ 




/ 


/ 

/ 


\\ > 


- / 


/ / 




\\\ 
\\' 








\\ \ 

v> 


r 






\\\ 

V\ 


- 






\\\ - 








\\ \ 

\\ - 








\\v - 








V\ 






1 1 





0.01 0.1 1 



k/k„ 

Fig. 12. — Spherically-integrated, compensated power spectra of 
the solenoidal component of velocity for driven strong-field MHD 
turbulence runs with fcpkL/27r = 2 at 10243 (solid), 512^ (long 
dashed), and 256^ (short dashed). The x-axis has been renormal- 
izcd to give fc/fcjv = 1- Also shown is a fit (dotted) to the slope of 
the 5123 run over the range 8 < fcL/27r < 18, P{k) oc fc-i-3B xhc 
10243 run appears only slightly steeper. 
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Fig. 13. — Spherically-integrated, compensated power spectra 
of the quantity u = p^/^v for driven strong-field MHD turbulence 
with fcpkL/27r = 2 at 10243 (solid). The X-axis has been renormal- 
ized to give fc/fejv = 1- Also shown is a fit (dotted) to the slope 
over the range 8 < kL/2iT < 36, P(fc) oc k'^ '^^. 

fcpki/27r — 8 shows more power at high wavenumber 
from Athena than from ZEUS for the velocity and its 
compressive and solenoidal components, just as was true 
for the hydro case. If the shallowing of our velocity spec- 
trum is due to the bottleneck effect, then what appeared 
to be inertial range when V03 was published is likely 
affected by the bottleneck as well. 

Figure shows our cylindrically-integrated power 
spectra at 512"^. We find anisotropy in the power spec- 
trum of magnetic field perturbations as well as the veloc- 
ity and its solenoidal component. There is more power 
perpendicular to the mean magnetic field at a given 
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Fig. 14. — Spherically-integrated, compensated power spectra of 
the compressive component of velocity for driven strong-field MHD 
turbulence runs with kpi^L/2n = 2 at 10243 (solid), 5123 (j^^g 
dashed), and 2563 (short dashed). The x-axis has been renormal- 
izcd to give k/kpf = 1. Also shown is a fit (dotted) to the slope of 
the 5123 run oygj. the range 10 < kL/2n < 18, P{k) oc fc-^ n. The 
10243 run, which does not have a sufficiently long time-average to 
be robust, does not appear to have power law form. The prominent 
upturn at the highest wavenumbers is due to the uncompensated 
power spectrum flattening out to a constant value. 
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Fig. 15. — Spherically-integrated, compensated power spectra 
of the magnetic field perturbations for driven strong-field MHD 
turbulence runs with kp]^L/2TT = 2 at 10243 (solid), 5123 (long 
dashed), and 2563 (short dashed). The x-axis has been renormal- 
ized to give k/kpf = 1. Also shown is a fit (dotted) to the slope of 
the 5123 run over the range 8 < kL/2iT < 18, P{k) oc fc-^'^^. 

wavenumber than parallel to it. Parallel to the mean 
field, the power law is roughly /c"^, similar to that of the 
purely hydrodynamic case. Perpendicular to the field, 
however, the slope is much more shallow, roughly 
We find the power spectrum of the compressive compo- 
nent of velocity to be nearly isotropic, in contradiction 
to what was found in V03. 

The velocity power spectrum of subsonic turbulence re- 
sults from a conservative cascade of energy from large to 
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Fig. 16. — Cylindrically-integrated (uncompensated) power 
spectra (solid) of (a) the velocity, (b) the magnetic field, and the 
(c) solenoidal and (d) compressive components of velocity for the 
driven strong-field MHD turbulence run at 512^ shown in Figure 
111! The X and y axes are k^^/kj^ and k±/ki^, respectively. The 
colors correspond to Px(fc|| , k±), where fc^ = fcy + k^. Also shown 
for reference are circular contours (dotted). 

small scales through interactions local in Fourier space. 
It is commonly believed that the velocity spectrum of su- 
personic turbulence, because it too is a power law, must 
also result from such a cascade. It is commonly known, 
however, that the power spectrum of a discontinuity, or 
an ensemble of shocks as in Burgers turbulence, also have 
velocity spectra that follow a power law, P{k) cx k^^. In 
Figure 1171 we confirm that the compressive component 
of velocity in our 256'^ supersonic, strong-field MHD tur- 
bulence run with fcpii;L/27r = 2 has a spectrum similar to 
that of velocity in an MHD blast wave that was initially 
spherical. This calls into question the long-held assump- 
tion that the power spectrum seen in supersonic turbu- 
lence results from an energy cascade facilitated by inter- 
actions local in Fourier space. It would seem that other 
diagnostics, such as structure functions, are necessary in 
order to determine the cause of the power law spectrum, 
either a Fourier-space cascade as in incompressible turbu- 
lence, or an ensemble of shocks as in Burgers turbulence. 

To get a clearer idea of which features in our power 
spectra are representative of the turbulence and which 
were introduced in the post-processing, we compare 
power spectra computed by multiple methods. Our con- 
trol spectrum will be computed as were the others pre- 
sented in this section — averaged over the cells within a 
shell and then multiplied by the volume of the shell, 
dV — Airk'^dk, which we will refer to as "shell-averaged" 
for the remainder of this section. Our comparison will be 
to another commonly used method (e.g. P07), where the 
power is simply summed over the cells falling within the 
shell ("shell-summed"). These methods differ due to the 
discretization of the grid, i.e. because the Fourier-space 
volume occupied by the cells within a bin is not, in gen- 
eral, equal to the volume of the perfectly spherical shell 
that the bin represents. While the bins in our control 
spectrum are centered on integer values of kL/2TT, those 
in the alternative spectrum instead run between integer 



0.1 - 




\ 



LiJ I I I I I I 

1 10 100 

k/k„ 

Fig. 17. — Spherically-integrated compensated power spectrum 
of the compressive component of velocity from one snapshot of 
the 256'^ driven strong-field MHD turbulence run (solid) shown in 
Figure 1111 compared to the total velocity power spectrum of an 
initially-spherical MHD blast wave (short dashed), also at 256^. 
Except for the oscillations, the shapes of these spectra look quite 
similar between the driving and dissipative scales. 

values. 

We use the time-averaged power spectrum of the driven 
strong-field MHD turbulence run at 512^ for our compar- 
ison, although the effect should be independent of the run 
being analyzed. When overplotting the spectra produced 
by these two methods (see Figure \TE\\ . the most obvious 
difference is at low k. Although the shell-averaged spec- 
trum is relatively smooth in the range 5 < kL/2'n < 20, 
the other has a jagged shape over this same range re- 
sulting from its sensitivity to the number of cells falling 
within a bin, making the "slope" of the power spectrum 
in this region much less obvious. 

If we compare the shell-summed spectrum to one com- 
puted by the same method but with the bins shifted by 
half a bin width (aligning these bins with our own) , we 
find that shape of the spectrum at low k changes consid- 
erably. If, on the other hand, we take the shell-averaged 
spectrum and compare it to one computed in the same 
manner but with bins shifted to align with those typ- 
ically used in the shell-summed method, we see that, 
although the power in each bin does change, the shape 
and slope of the spectrum change very little. While nei- 
ther method is right or wrong, we advise caution when 
choosing a binning method for power spectra; differences 
in slope computed by different means are not necessarily 
indicative of different turbulent states. 

3.3. Sonic Scale 

Although molecular clouds have supersonic velocity 
dispersions on large scales, as one looks to smaller and 
smaller scales, turbulent compressions will become pro- 
gressively weaker, at some point becoming subsonic. The 
length scale at which the RMS velocity dispersion is equal 
to the sound speed is referred to as the sonic scale (Mc- 
Kce & Ostriker 2007). We now investigate how the veloc- 
ity dispersion varies on spatial scales larger than the sonic 
scale, where shocks are most important. This scaling can 
be determined observationally in the form of linewidth- 
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Fig. 18. — Spherically-integrated compensated power spectrum 
of the total velocity for t he 512^ driven strong-field MHD turbu- 
lence run shown in Figure lTTI calculated using our standard method 
(solid) compared to the method described in P07 (dotted). The 
PS computed using the P07 method is more jagged in the range 
5 < kL/2TT < 20. Note that the compensation used for this pair of 
power spectra is non-standard. 



Fig. 19. — 3D velocity dispersion versus the spatial scale on 
which it is measured for driven hydrodynamic turbulence with 
fcpijL/27r = 2 (open hexagons). Also shown is a power law fit, 
crl oc fi-^^, from the driving cutoff {kL/2iT = 8) down to our last 
data point above the sonic scale (both limits marked with dashed 
lines) . We have connected our data points with a solid line to make 
comparing to the fit easier. 



size relations using many different methods. For ex- 
ample, Brunt (2003) used principal component analysis 
(PCA) to determine the linewidth-size relation within 
individual clouds. 

To determine the velocity dispersion at a given length 
scale, fcL/27r = m, where m = 2" and n is an integer, 
we divide our computational domain along each axis into 
m sections, yielding sub-volumes. We compute the 
velocity dispersion in each sub-volume and then aver- 
age over all sub- volumes at that scale. Because crossing 
the sonic scale represents a change in the physics dom- 
inating the flow, i.e. waves steeping to form shocks, for 
this analysis only we will compute the velocity dispersion 
without mass- weighting, i.e. — [a^^ + +o'^J^^^i 

where a^^ = [{vf) - {vi)"^]^/^. 

Figure [19] shows the velocity dispersion versus length 
scale for the 1024^, M 7 driven hydrodynamic turbu- 
lence run with fcpki/27r — 2. Since driving may affect 
the scaling relation, we truncate our driving spectrum 
at kL/2'K = 8, where the driving has already dropped 
to only 1% of the peak level, and consider only this and 
smaller scales. Fitting to the data points falling between 
the sonic scale and the driving cutoff, a factor of eight in 
length, we find a clear power law of the form v{l) oc f'-^^ . 
This scaling index falls well within the range of indices 
found observationally by Brunt (2003). 

Shown in Figure [20] is the velocity dispersion versus 
length scale for the 1024^, M. ^ 7 driven strong-field 
MHD turbulence run with kpiiL/2'K ~ 2. This time, over 
the range of scales where the flow is supersonic, we find 
that the slope is not a power law. Unlike in hydrody- 
namic turbulence, the sonic scale is not the only scale at 
which the dominant physics changes. In fact, we would 
expect the multiple wave families of MHD to lead to mul- 
tiple transitions, and for strong fields, these transitions 
are widely separated. For example, parallel to the mag- 
netic field, slow waves will travel at the sound speed. 




0.01 0.1 1 

Length scale 



Fig. 20. — ID velocity dispersion parallel to the magnetic field 
(open squares) and 2D velocity dispersion perpendicular to the field 
(open triangles) versus the spatial scale on which it is measured for 
driven strong-field MHD turbulence with fcp]jL/27r = 2. From the 
driving cutoff {kL/2-K = 8) down to our last data point before the 
3D velocity dispersion becomes subsonic (both marked with dashed 
lines), neither the parallel nor perpendicular velocity dispersion has 
power law form (dotted). We have connected our data points with 
a solid line to make comparing to the power laws easier. 

whereas fast waves will travel at the Alfven speed. Per- 
pendicular to the field, however, the slow waves will have 
zero velocity and the fast waves will travel with a speed 
\/v\ + cl- This is complicated further by the fact that 
the parallel and perpendicular directions are defined rel- 
ative to the local magnetic field, not the mean field. As 
a result, we have no reason to expect a uniform power 
law between the sonic scale and driving cutoff. It is not 
clear to what extent this result is affected by dissipation. 
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3.4. Time-variability and Temperature Intermittency 

Up to this point we have averaged quantities from mul- 
tiple snapshots in order to minimize the effect of rare 
events on our results. Our final set of diagnostics will 
instead be an analysis of the time-variability and tem- 
perature intermittency of the turbulent gas. 

3.4.1. Saturation Amplitude 

We begin by analyzing the time-variability of some of 
the quantities we analyzed in ^3.11 For our hydrody- 
namic runs with the smallest driving scale, fcpki/27r = 8, 
we find the time-variability of the saturation energy to 
be only half a percent at both 256^ and 512^. Time- 
variability of less than a percent agrees with the results 
in S98 even though they drove their turbulence impul- 
sively. The time-variability quoted in K07, however, is 
much higher. The K07 runs were driven at a much larger 
driving scale than our own. To investigate if this could 
be the reason for the discrepancy, we will also analyze the 
time- variability of our runs with larger driving scales. 

When we increase the driving scale to fcpki/27r = 4, 
we find the time- variability of the energy to be a full per- 
cent at both 256'^ and 512^. As both the driving scale and 
Mach number are varied here as compared to the run we 
discussed above, this could be due to either. When we in- 
crease the driving scale further, to fcpki/27r — 2, without 
changing the Mach number, we find the time-variability 
increases yet again, coming to 2% at 256^ and approach- 
ing 4% at 512'^. It would seem that the driving scale 
chosen has a significant impact on the time- variability of 
the energy. The values we find here are still less than 
that shown in the K07 figure, however. 

For strong-field MHD, we again find time- variability 
of about half a percent for the fcpki/27r = 8 runs. In- 
creasing the driving scale to kp^L /2Tr — 4 decreases 
the time-variability at 256'^ but increases it a small 
amount at 512^. Increasing the driving scale further, to 
A;pkl'/27r = 2, yields time-variability of 3% at 256^ and 
4% at 512^, comparable to the hydro runs. 

When we compare the time-variability of the energy 
among hydro runs at 512"^ with varied Mach number that 
are driven at fcpk-^/27r — 4, we find that the relation 
is non-monotonic. With fcpki/27r — 2, however, we do 
see a monotonic relationship, with the time-variability 
increasing with Mach number. For our strong-field MHD 
runs, we again find non-monotonicity at fcpki/27r = 4, 
while the time-variability increases with Mach number 
at fcpki/27r = 2, with the exception of the lowest Mach 
number run. 

3.4.2. Temperature Intermittency 

A feature of turbulence that may have a strong in- 
fluence on star formation is intermittency, dissipation 
that is highly localized in space and time (McKee & 
Ostriker 2007). Although such dissipation does not re- 
quire shocks, it is worthwhile to study the shock con- 
tribution in supersonic turbulence. Were we not assum- 
ing an isothermal equation of state, energy dissipation 
in these regions would lead to local heating and thus in- 
creased temperatures, which would be evident in observa- 
tions. Therefore, to study intermittency due to shocks, 
we analyze the maximum heating rate per unit mass, 
Q — c^V • V, a surrogate for temperature. 



Since we expect the maximum heating rate, Qmax, to 
be strongly influenced by rare events, such as the interac- 
tion of multiple shocks, we consider the time-variability 
of this quantity. Due to the discretization of the grid, 
however, we would expect a component of the time- 
variability to be due to grid-scale fluctuations. There- 
fore, instead of considering Qmax directly, we consider 
the statistics of the high-Q tail of the PDF of the heat- 
ing rate per unit mass. For each simulation that we con- 
sider, we compute the heating rate per unit mass, Q, 
in each cell over many snapshots. Considering only the 
cells in each snapshot that compose the 1% of material 
(by mass or volume) with the highest cell-averaged Q, we 
compute the mean value. Following the time-variability 
of this value should allow us to ignore meaningless grid- 
scale fluctuations while studying the intermittency due 
to shocks. 

For driven hydrodynamic turbulence at 512^, we find 
the low Mach number runs driven at k-^yh /2ti = 4 to have 
less than 1% variability in the peak temperature, while 
the higher Mach number runs with this same driving 
scale can have variability as large as 2% but without a 
monotonic dependence on Mach number. At the larger 
driving scale fcpki/27r = 2, however, the variability in the 
peak temperature always exceeds 1% and can be as high 
as 4%. It would seem that this, like the time-variability 
of the saturation energies, increases with driving scale. 
For our runs, the variability of the peak temperature as 
measured by the top 1% of the volume never exceeds that 
from the top 1% of the mass by more than a tiny amount. 
As we have very few data dumps with which to calculate 
the statistics, however, the quantitative behavior may 
not be robust. 

We find that, in 512^ strong-field MHD turbulence 
driven at fcpki/27r = 2, the variability of the peak 
temperature can be as large as 9%. When driven at 
k-p\^L/2'K = 4, however, the variability always decreases, 
never exceeding 3%. The variability measured by the 
top 1% of mass is typically larger than that from the 
top 1% of volume, but there does not appear to be a 
strong dependence on Mach number. More often than 
not, the variability in the strong-field MHD runs are 
larger than in the hydro runs at the same Mach number. 
Again, however, these statistics are computed from very 
few data dumps, making them subject to large errors. 
Figure [2T] shows the time-evolution of the peak temper- 
ature in the fcpki/27r = 4 strong-field MHD run with 

E/pL'^Cg = 1000, for which we have better than typi- 
cal statistics, compared to the fcpki/27r = 2 run with 
E/pL^cl = 500. 

4. DISCUSSION AND CONCLUSIONS 

The saturation energies and dissipation timescales we 
find support the conclusion of S98 that supersonic tur- 
bulence dissipates rapidly even in the presence of mag- 
netic fields. The results of our Godunov code agree with 
ZEUS, an operator-split method that relies on artificial 
viscosity to capture shocks, but reach convergence at 
slightly lower resolutions. At 512^ the difference in the 
total energy in fluctuations between Athena and ZEUS 
for strong-field MHD is very small, indicating that these 
two codes converge to similar turbulent states at suffi- 
ciently high resolution. 

The convergence resolution of our simulations depends 
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Fig. 21. — Maxi mum heating rate per unit mass (solid), where 
the tail (see ^3.4.21 1 is defined hy the top 1% of mass, for driven 
strong-field M ~ 7 MHD turbulence with fcpkL/27r = 4 at 512^. 
The time-variability here is lower compared to the case where 
kp]^L/2iT = 2 (shown by stars) at roughly the same Mach num- 
ber. 

both on the driving scale and the presence of a magnetic 
field. While hydro turbulence driven at fcpki/27r = 4 
converges by 64^ , strong- field MHD turbulence driven at 
fcpki/27r = 8 does not converge until 512^. Although 
very high resolutions are needed to capture the inertial 
range of the turbulence, lower resolutions are often ad- 
equate for studying the energy dissipation characteris- 
tics of the turbulence as well as other volume-integrated 
quantities. 

At high Mach number, the ratio of the dissipation 
timescale to the flow crossing time at the driving scale in- 
creases with increasing magnetic field strength; however, 
it does not exceed unity even for strong-field MHD tur- 
bulence. This ratio is independent of driving scale. The 
fractions of the kinetic energy in solenoidal and compres- 
sive modes are also independent of driving scale, with the 
compressive fraction generally being more than twice as 
high for hydrodynamic than for strong-field MHD turbu- 
lence. 

Our spherically-integrated velocity power spectra for 
decaying, subsonic hydrodynamic turbulence show evi- 
dence of the bottleneck effect in the velocity power spec- 
trum, consistent with the findings of Sytine et al. (2000). 
We find more power at high wavenumber in our driven, 
supersonic hydro and MHD turbulence simulations than 
was found by V03, but it is unclear, particularly in the 
MHD case, whether or not this is due to a bottleneck. 
Resolutions exceeding 1024^ will be necessary to draw 
firm conclusions about the slope of the inertial range. 
The cylindrically-averaged velocity power spectrum for 
driven MHD turbulence is very anisotropic; it has a slope 
that approximates that of the hydrodynamic case paral- 
lel to the magnetic field, while perpendicular to the field 
it is much more shallow. The compressive component 
of velocity has an isotropic power spectrum, contrary to 
what was found in V03. 

We find the compressive component of velocity in 
driven, supersonic MHD turbulence to have a power 
spectrum that is difficult to distinguish from the veloc- 



ity spectrum of an initially-spherical MHD blast wave. 
This calls into question the long-held assumption that 
supersonic turbulence power spectra result from an en- 
ergy cascade facilitated by interactions local in Fourier 
space. The analysis of structure functions may be useful 
in determining the source of the power law spectrum in 
supersonic turbulence, either a Fourier-space cascade as 
in incompressible turbulence, or an ensemble of shocks 
as in Burgers turbulence. 

For hydrodynamic turbulence, we find a power law 
scaling of the velocity dispersion with spatial scale, 
a{l) DC 1°-^^, for scales where the velocity dispersion is su- 
personic. However, we find no such power law for strong- 
field MHD turbulence, where the sonic scale is not the 
only scale of interest. In this case, we find that the ve- 
locity dispersion drops off more rapidly than a power law 
as one approaches smaller scales. 

We see time- variability in the saturation energies com- 
parable to that of the equivalent runs in S98 despite the 
impulsive driving employed therein. Our time-variability 
increases when we apply our turbulent driving at larger 
scales, but remains lower than that shown in K07 even 
with fcpkL/27r = 2. It is possible that the acceleration- 
based driving method of K07 is responsible for the dif- 
ference. At this large driving scale, the time-variability 
increases with Mach number for both hydro and strong- 
field MHD. The method used to drive the turbulence 
appears to have a substantial impact on the resulting 
turbulent state. 

Further investigation should be conducted to deter- 
mine how much of the increase in time-variability with 
driving scale is due to the limited range of scales over 
which an inverse cascade can occur. If the level of time- 
variability is determined to be a result of simulation 
setup and not representative of the physical system we 
are trying to simulate, better diagnostics or simulation 
methods need to be developed in order to quantify in- 
termittency in turbulent media. Although we have very 
poor statistics in our temperature intermittency anal- 
ysis, it appears that this, too, increases with driving 
scale. Considering the difference that the driving method 
makes on the results, it may be more realistic to study 
decaying turbulence instead of driving it arbitrarily. 

For identical turbulent data cubes, we find that the 
post-processing methods chosen significantly influence 
the results. The turbulent Mach number changes by 
~ 4% depending on the method used in its computa- 
tion. The means of computing the power spectrum also 
has a strong influence on the power at wavenumbers over- 
lapping with the inertial range. When comparing results 
published by different groups, one should keep in mind 
that this, as well as the code, driving method, and initial 
conditions, can affect the quantities being compared. 

Although the saturation energies and energy dissipa- 
tion characteristics of the turbulence converge at resolu- 
tions within our current computational capabilities, the 
power spectra appear to require much higher resolutions 
to provide valuable information. Also considering that 
turbulent power spectra can be approximated by non- 
turbulent phenomena, it would seem that, for the time 
being, our focus should be on other diagnostics. 

In conducting these numerical simulations, we have 
made many simplifications in order to make the problem 
more tractable. These assumptions, however, may prove 
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to significantly impact the results. Future studies should 
consider non-ideal MHD in order to model low-ionization 
molecular clouds, a non-isothermal equation of state in 
order to study heating and cooling, and self-gravity in 
order to follow the collapse of the bound clumps that 
form in the turbulent medium. 

We thank Eve Ostriker for very productive discussion 
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Kritsuk, Cristoph Federrath, and Alex Lazarian for their 
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Gene at Princeton and on computational facilities sup- 
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APPENDIX 

SUPERSONIC TURBULENCE WITH GODUNOV SCHEMES 

The calculations presented in this paper were conducted using Athena, a directionally-unsplit, higher-order Godunov 
code. This code exactly conserves mass, momentum, and magnetic flux, as well as energy when assuming an adiabatic 
equation of state. The code captures shocks well and has a low level of numerical dissipation. Although the full details 
of the algorithms and implementation can be foimd in the literature (Gardiner & Stone 2005, 2008; Stone et al. 2008; 
Stone & Gardiner 2008), we will briefly summarize here, noting any modifications made to the algorithms in order to 
successfully run the challenging problem of high Mach number turbulence. 

The integration algorithm used in our calculations is a simple predictor-corrector scheme based on the MUSCL- 
Hancock scheme of van Leer (2006), combined with the constrained transport method of Evans & Hawley (1998) 
to enforce the divergence-free constraint on the magnetic field (i.e. the VL-I-CT algorithm described in detail in 
Stone & Gardiner 2008). We find the additional diffusion associated with this scheme as compared to our CTU+CT 
algorithm (described in detail in Stone et al. 2008) to make it more robust to the strong rarefactions that arise 
in a highly turbulent medium. Although the algorithm is formally only second-order accurate, we use third order 
(piecewise parabolic) spatial reconstruction, finding it to provide more accurate solutions in test problems due to 
smaller truncation error. 

For our isothermal hydrodynamics runs, we found there to be strong rarefactions within the turbulent medium for 
which an approximate Riemann solver simply was not accurate enough, necessitating the use of an exact nonlinear 
Riemann solver. For adiabatic hydrodynamics as well as isothermal and adiabatic MHD, we were able to use approx- 
imate nonlinear Riemann solvers, namely HLLC for the hydro case and HLLD for MHD, more details about which 
can be found in Stone et al. (2008). Although we used our own adaptation of the adiabatic HLLD solver of Miyoshi & 
Kusano (2005) for our isothermal MHD runs, we found it to produce turbulent states extremely similar to those from 
the isothermal HLLD solver of Mignone (2007). 

Although these Riemann solvers are positive definite in ID, it is not guaranteed that they will be so in multidimen- 
sions. In fact, we found that, under extreme conditions, the net mass flux out of a cell in our isothermal 3D turbulent 
medium occasionally exceeded the cell's initial mass. In the rare circumstance that this occurred, we recomputed 
the fluxes of all variables through all interfaces bordering such cells using first-order reconstruction. In our Mach 7 
strong-field MHD run at 512"^ that uses fcpkL/27r = 2, this affected only a fraction 3 x 10^^" of the fiuxes computed in 
the corrector step of the integrator. Dropping to first order introduced enough diffusion in the immediate vicinity to 
keep the cell-averaged density positive, while having a negligible effect on the overall system. Adding diffusion in this 
manner instead of enforcing a density floor maintains exact conservation. 
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